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Abstract: The goal of Point Distance Solving Problems is to find 2D or 3D placements of points 
knowing distances between some pairs of points. The common guideline is to solve them by a 
numerical iterative method {e.g. Newton-Raphson method). A sole solution is obtained whereas 
many exist. However the number of solutions can be exponential and methods should provide solu¬ 
tions close to a sketch drawn by the user. Geometric reasoning can help to simplify the underlying 
system of equations by changing a few equations and triangularizing it. This triangularization is 
a geometric construction of solutions, called construction plan. We aim at finding several solu¬ 
tions close to the sketch on a one-dimensional path defined by a global parameter-homotopy using 
a construction plan. Some numerical instabilities may be encountered due to specific geometric 
configurations. We address this problem by changing on-the-fly the construction plan. Numerical 
results show that this hybrid method is efficient and robust. 

Key-words: Point Distance Solving Problems, Reparameterization, Curve Tracking, Symbolic- 
Numeric Algorithm 
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Une methode robuste et efRcace pour resoudre des 
systemes de contraintes de distances entre points par 

homotopie 

Resume : Le but de la resolution de problemes de contraintes de distances entre points est de 
placer en 2D ou 3D un ensemble de points connaissant certaines distaces entre paires de points. 
De tels problemes sont en general resolu grace a une methode numerique, souvent Newton- 
Raphson, qui ne produit qu’une solution alors qu’il en existe un nombre exponentiel. Cedes 
ressemblant a I’esquisse sont d’un interet particulier. Le raisonnement geometrique pent cepen- 
dant aider a simplifier les sytemes d’equations correspondant en remplagant quelques equations, 
ce qui permet de les triangulariser. Une telle triangularisation est une construction geometrique 
des solutions et est appellee plan de construction. On se propose dans ce rapport de trouver 
plusieurs solutions, proches de I’esquisse sur une courbe definie par une homotopie utilisant le 
plan de construction pour reduire son cout. L’utilisation d’un plan de construction induit des 
instabilites numeriques a proximite de certains points; ces instabilites sont evites en changeant 
le plan de construction pendant le suivi de la courbe. La methode decrites ici a ete implementee, 
et les resultats obtenus montrent son efiicacite et sa robustesse. 

Mots-cles : Problemes de Constraintes de Distances entre Points, Re-parametrisation, Suivi 
de courbes, Algorithme symbolique-numerique 
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1 Introduction 

Geometric Constraints Solving Problems arise in many fields such as CAD, robotics or molecular 
modeling. The problem is to determine the positions of geometric elements (points, lines, planes, 
circles, etc.) that must satisfy a set of constraints such as distances, angles, tangencies and so 
on. Commercial solvers generally rely on numerical methods such as Newton or quasi-Newton 
that provide a single solution. Even when problems are well-constrained the number of solutions 
may grow exponentially with the number of constraints. Usually the user is not interested in all 
solutions but only to those whose shape is close to a sketch that he provided. 

Here the restrained class of problems involving only points and distance constraints is consid¬ 
ered. Our aim is precisely to design a method that uses the sketch to guide the research of several 
solutions. We assume that the considered problems are structurally well-constrained. Roughly 
speaking, this means that there exist some assignments for dimensions leading to hnitely many 
solutions. We also consider problems resisting to decomposition-recombination methods. 

Several methods can yield several or even all the solutions. Subdivision methods [m E] 
provide all the solutions but the number of boxes to be explored can be huge. In algebraic 
approaches, homotopy methods have been successfully studied in this area lain] but only for 
small size problems. Indeed, the number of homotopy paths to follow grows exponentially with 
the number of constraints. In m it is proposed to use the sketch to define a parameter-homotopy. 
A sole path is followed but a sole solution is obtained. 

Another way to get several solutions comes from geometric methods. A construction plan 
is first derived by applying some geometric construction rules. It consists in a sequence of 
basic construction steps. Next, such a plan is numerically evaluated to yield different solutions. 
However, no construction plan can be easily found for some 2D problems and for most 3D 
problems. To circumvent this, proposes a new approach that performs a reparameterization. 
In this approach, a geometric constraint system, say ^i, is modified by adding and removing 
some constraints to obtain a system similar to the original one and from which a construction 
plan can be easily derived. In turn, this construction plan is used to define a reduced system R 
from the constraints removed from Si. i? is then solved by a numerical solver in order to meet 
the removed constraints while still satisfying the constraints of S 2 . For point distance problems, 
the number of constraints that have to be swapped to obtain S 2 from Si is much lower than the 
number of constraints of hence the size of the system to be numerically solved is drastically 
reduced. 

The drawback is that the equations of system R are much harder to deal with due to irregular 
configurations. So the choice of the numerical method is crucial to provide several solutions. In 
[7] one or two constraints are removed and the solutions are found by a sampling method. In [3] , 
this idea is extended for more than two constraints, Newton-Raphson method allows to get some 
of the sought solutions. In [5], the reparameterization is used at a low-level to simplify linear 
algebra involved by numerical methods. Finally [5] presents a first attempt for using a homotopy 
method along with reparameterization. The idea is to follow a homotopy path to which belongs 
the sketch. 

All these work can quickly find some of the solutions desired by the user. However some 
solutions are often missed because of numerical inaccuracies. In this paper we provide an effective 
and original method to face it. This work is based on tracking homotopy paths defined by a 
construction plan obtained after reparameterization of the problem (system S 2 ). The central 
idea is to detect ill-conditioned configurations induced by the interpretation of the construction 
plan, and to change on-the-fly the construction plan to get away from such configurations. We 
justify this approach by showing that these changes of construction plans during paths tracking 
do not change the path that is followed. More precisely. 
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• we show that the paths followed with a homotopy method applied to the reduced system 
using a construction plan can be glued together around the singular points. The whole 
path is exactly the path which would be followed by continuation on the original system. 
Thus, it is independent of the reconstruction plan obtained with the reparameterization 
phase. 

• we use geometric criteria to detect in advance the singular points caused by a particular 
construction plan, and we design a way to modify on-the-fly the reparameterization in 
order to avoid singular points. 

• putting everything together, we marry a homotopy method with a reparameterization to 
have a new algorithm to solve point distance satisfaction problems. We prove that this 
algorithm terminates and is correct. We compare the results obtained with our new method 
with another homotopy method and we find that this is algorithm is more than three time 
faster than our previous algorithm without reparameterization. 

The rest of the paper is organized as follows. Focusing on 2D problems, Sec.j^gives definitions 
on construction plans and homotopy. Sec. gives results about homotopy paths tracking on 
construction plans that justify our approach. Sec. [^explains how to change a construction plan 
on-the-fly to avoid critical situations. The soundness of the approach is justified in Sec. Sec. 
gives tracks to extend our method to 3D problems, and Sec. [^presents some experimental results. 


2 Notations and Definitions 

A Point Distance Satisfaction Problem (PDSP) is a constraint satisfaction problem where con¬ 
straints are imposed distances between points. Unknown points are sought either in the Euclidean 
plane for 2D PDSP or in the Euclidean space for 3D PDSP. We focus here on the 2D case. 

The method presented in this paper uses symbolic manipulations on PDSP to ease their 
numerical solving. For the sake of clarity in the description of this symbolic-numeric approach, 
we will use different typefaces to denote a variable and its numeric value. A boldface lowercase 
letter as x will denote a variable, and an uppercase boldface letter as X will denote a set of 
variables. A value for x, in general a real number, will be denoted by a lowercase italic letter x, 
and a value for X, in general a real vector, will be denoted by a uppercase italic letter X. We 
make an exception for variables associated with geometric objects: if p is a point, we will note p 
a value for p whereas it refers to a vector of real values (its coordinates) in a geometric context. 

2.1 Point Distance Satisfaction Problems 

A PDSP G is denoted by G = G[P, A] where P is a set of unknown points, A is a set of length 
parameters, and G is a set of m constraints of distance. A distance constraint of parameter 
ai S A between points P2,P3 G P is written distance{p 2 ,P 3 ) = ai. A PDSP that consists in 
constructing 6 points pi, ... ,P 6 in a plane knowing 9 distances is given in Fig. We call it 

A 3 Numerical values for A are usually given by a user. The aim is to find the solutions 
that respect these dimensions. We suppose in addition that a sketch, i.e. a geometric placement 
of points of P, with possibly a representation of constraints, is available. Right part of Fig. 
shows a sketch of the PDSP ^ 3 ^ 3 . 

^when considering right part of Fig.j^as a non-oriented graph, it is the complete bipartite graph with 3 vertices 
in each component 
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Unknowns: 

point pi,pe 
Parameters: 
length ai,ag 
Constraints: 

distance(pi, P2) = ai 
distance{p 2 ,P 3 ) = a.2 

Pi ai P2 

distance{p2,ps) = as 
distance{p3,ps) = ag 

Figure 1: A symbolic statement (left part) and a dimensioned sketch (right part) of the PDSP 
A' 3 ^ 3 . Edges are distance constraints of parameters a^. 



Unknown points Pj are sought in the Euclidean plane and are each associated with two 
algebraic unknowns corresponding to their coordinates. Let Ci be distance{pj,pk) = a,;. It is 
associated with the numerical function Ci(P, A) called numerical interpretation: 


Cj(P, A) = PjPfe - ai (1) 

where PjPk holds for the Euclidean distance between pj and p^. Since constraints of distance are 
invariant up to rigid motions of the plane, placements of points in the plane fulfilling constraints 
are sought in a reference, i.e. the values of 3 unknown coordinates are fixed. For we could 
search values for points with pi at the origin and P 2 with null ordinate, and assign variables 
{xi, ...,Xg} to remaining free coordinates. 

We denote by X = {xi, ...,Xm} the set of free unknown coordinates of points in P and we 
define the system of equations associated to G as 


F(X,A) =0 


(■^) 


where F : K™ x —>■ K™ has as i-th component the numerical interpretation of the constraint 

Ci defined in Eq. 0. 

We will call figure a set of real values X for X. Given positive values for A, we call 
solution of G a figure X that is a solution of [J^sol defined as 

E(X,A«°)=0 (J-^o) 


We highlight here that a sketch of G is a figure In addition, by measuring on A®* 

distances between appropriated points one can find positive values A®^ for A such that A®^ is a 
solution of the system Esk defined as 


F(X,A®'=) = 0 (Esk) 

In the following we will consider structurally and generically well constrained PDSP. A struc¬ 
turally well constrained PDSP satisfies in particular |X| = |A| = |G| = to if elements of A 
are algebraically independent (see the Koenig-Hall theorem [H]). A generically well constrained 
PDSP admits for generic values of parameters a not null and finite number of solutions. Here 
generic stands for the complementary of a set having a null Lebesgue measure in an open subset 
of the space of parameters. 
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2.2 Homotopy 


Equations of |J^so| can be written as polynomials, hence [7^ can be solved in C”^ by a classical 
homotopy method (see [T] for an introduction to homotopy methods, and [3] for its application 
to our context). All the complex roots ofare searched and found, whereas in applications 
only the real solutions are relevant. Here we aim at obtaining only real solutions ofand we 
use the sketch to define a real homotopy (see [11113]) between Tsk and Pgo using an interpolation 
of parameters defined as follows. 


Definition 1 Let = {af^, ..., a^} and A®° = {af°,..., a™} be strictly positive real values 
for A. We call interpolation function from to A®° a C°° function a : K —>■ K"* that satisfies: 

• for 1 < i < m, ai{0) = af^ and 0 ^( 1 ) = a®°, and 


• for 1 < z < m, Vt e [0, l],ai{t) > 0. 

We call positive support of an interpolation function a and we note it supp+(a) the subset o/K 
where Oift) > 0 for 1 < i < m. 

Notice that [0,1] C supp+{a). Given an interpolation function a from A®^ to A®°, we define 
the homotopy system jW] as: 

iJ(X,t)=F(X,a(t)) = 0 {%) 

where PI : M*” x K —)■ K™ is called the homotopy function associated to G. 

The set of solutions of \m\ denoted by H“^(0) can be partitioned into a set of connected 
components that are called homotopy paths of It is worth mentioning here that a point {X, t) 
with real components belongs to a homotopy path of only if t S supp+{a), otherwise two 
points of X would be separated by a negative length. The following result (see |9]) underlines 
the influence of a and supp+{a) on the topology of homotopy paths of|^ It is here stated in the 
general case where constraints are not necessarily distances but also angles, collinearities and so 
on. 


Theorem 1 ([9]) Let H define the homotopy^^ Jh be its Jacobian matrix, PH C x K 6e 
its domain of definition, and PH‘^ be its complementary. Under assumptions 

(hO) PH is open, 

(hi) Jh has full rank on each point of H~^{0), 

(h2) supp+{a) is compact, 

homotopy paths o/H are 1-dimensional manifolds diffeomorphic either to circles, or to open 
intervals. If an homotopy path S is diffeomorphic to an open interval then the extremities of S 
converge either to a point in PH'^ or to a solution with an infinite norm. 


We make here two remarks to adapt this result in the framework of PDSP. First, H is clearly 
C°° in X K, hence paths cannot converge to a point of PH^. Secondly, if supp+{a) is compact, 
a{supp+{a)) is compact and all components of a(t) are bounded if t S supp+{a). Hence if {X,t) 
belongs to an homotopy path ofj^ the components of X are coordinates of a set of points lying 
in a compact. In a PDSP context, Thm.l^can be restated as follows: 

Corollary 1 Let G be a PDSP and H define the homotopy satisfying assumptions (hi) and 
(h2) of Thm. Homotopy paths of^ are 1-dimensional manifolds diffeomorphic to circles. 

Denoting by S the homotopy path to which belongs (A®^,0), S can be followed with a 
numerical path tracker until it loops on (A®*, 0). It allows to find points {X,t) of S with t = 1 
that are real solutions of [Asp [ 
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Unknowns: 

point pi,pe 
Parameters: 
length ai,ag, k 
Constraints: 

distance(pi, P 2 ) = ai 
distance{p 2 ,P 3 ) = a.2 
distance{pi,p3) — k 

distance{p2,P3) = ag 

Figure 2: A symbolic statement (left) and a dimensioned sketch (right) of the PDSP qK^ ^. 

2.3 Reparameterization 

A PDSP can be solved very easily when its associated system|^can be organized in a triangular 
form. From a geometric point of view, solving such a system is done by constructing points 
iteratively as intersections of two circles (three spheres in a 3D context) while making choices 
between possible intersections. The formal statement of the latter geometric construction is 
called a construction plan. When a PDSP G cannot be solved with this approach, an idea 
called reparameterization (see H) is to introduce d new constraints called added constraints 
with unknown parameters called driving parameters, in such a way that a construction plan 
parameterized by driving parameters constructs hgures fulfilling all constraints but d that are 
called removed constraints. G is then solved by hnding values of driving parameters such that 
the figures constructed by the construction plan satisfy the removed constraints. 

2.3.1 Construction Plans 

Consider the PDSP depicted in Fig. {q holds for “quasi”) that has been obtained from 

^ 3,3 by substituting the constraint distance{p^,pQ) = ag by distance[pi,p^) = k. Knowing 
values for {ai,..., ag, k}, its solutions are all found by the simple ruler and compass construction 
given in the leftmost part of Fig. 

The instruction pg = /nter(7L(pi, ai, li) holds for the construction of P 2 as one of the 
intersections of the line li with a circle of center pi and radius ai. Here pi and li are objects 
of the reference that are fixed to construct solutions up to rigid motions. The instruction pi = 
InterCG{pii , aj 2 , pp , ) holds for the construction of Pi as one of the intersections of the circles 

respectively centered in pp and pp of radius ap and ap. 

We will note /i[pi+i, A^] the instruction that constructs Pi+i from objects Ai (after a pos¬ 
sible re-indexing of points). Notice that objects of Ai are not only length parameters but also 
geometric objects of the reference or objects constructed by previous instructions. 

Definition 2 (CP) A Construction Plan (CP) of objects P and parameters A with reference 
Aq is a finite sequence I = of terms /i[pi+i,Ai] in a triangular form, i.e. 

(i) each p € P is either in Ag or is constructed by an instruction C S I, 

(a) for each instruction Ajpi+i, A^] G I, each a G Ai is either in Aq, or in A, or is constructed 
by a term Ij G I with j < i. 

We note it /[P,A,Ao], or more simply I. 

A CP can be seen as a symbolic solution of a set of constraints. Consider for instance an 
instruction Pi+i = /nterC'C'(pp, ap , pp, ap), it gives a symbolic solution to the constraints 
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Unknowns: 

point p 2 ,pe 

Parameters: 

point pi, line li 
length ai,ag, k 

Terms: 

P 2 = /nterC'L(pi, ai, li) 

P3 = /nterC'C'(pi,k, p 2 ,a 2 ) 
P 4 = 7nterC'C'(pi,a7,p3,a3) 
ps = 7nterCC'(p2,a8,P4,a4) 
P6 = 7nterC'C'(p5,a5,pi,a6) 



Figure 3: Left: a construction plan of solutions of the PDSP qK^ ^. Middle: a sub-tree of the 
associated interpretation tree. Right: two figures obtained when evaluating the CP on the two 
branches in solid lines. 


distance(jpi+i,Pii) = 3.^2 and distance{pi+i,Pi3) = We associate in such a way a set Cj of 
constraints with a CP I, and we say that 7[P, A, Aq] is a CP of a PDSP G = Cp, A] if Cj = C. 

2.3.2 Evaluation of a Construction Plan 

Given values A and Aq for A and Aq, a CP 7[P, A, Aq] is evaluated to obtain numerical values 
of the solutions of constraints C/p, A] by sequentially applying its instructions. At each step 
a choice between two intersections is done and considering all the possible intersections leads to 
construct an interpretation tree; its branches bring numerical values for P. Middle part of Fig. 

shows a sub-tree of the interpretation tree associated to the CP presented in the left part, and 
in its rightmost part it shows the two figures brought by the two branches in solid line. 

Let 7i[pi_|_i, Ai] be an instruction of 7. It is interpreted by a multi-function that maps to a 
value Ai of Ai the two possible intersection locii of objects defined by Ai. We index these locii 
by an integer and for a given index bi we note ^^'^Ii{Ai) the function that maps to Ai the locus 
of index bi. We consider that is not defined when the number of intersections is zero or 
infinite. We assume that indexation of intersections is continuous, i.e for each 1 < i < I and for 
each index bi, ^^'^Ii{Ai) is C°° on the interior of its domain of definition. 

We call branch of 7 a sequence b = {bi)\^i of indexes and, for a given branch b, we call 
evaluation of I on its branch b the numerical function A) that maps to values {Aq,A) 

the composition of functions On a given branch b, ^^^7 is C°° in the interior of its domain 

of definition as a combination of (7°° functions. 

2.3.3 Reparameterized Construction Plans 

Reparameterized construction plans (RCP) are central objects in the method presented in this 
paper. They appear in two steps. First, each problem must be derived in a CP. We are interested 
in this article in problems whose construction is not known, it is then necessary to transform the 
problem by adding and removing constraints as explained above. Constraints are added so that 
a CP is easy to establish. We do not detail here the way a RCP with a sole driving parameter in 
each instruction is obtained, see for different approaches. This new CP is called a RCP 

and is completely characterized by: the CP itself, the removed constraints that must be satisfied 
by all solutions and the driving parameters that are the added dimensions. For instance, a RCP 
for could be given by the CP 7 given in fig. the removed constraint (iis7(p3,p6) = ag 
and the driving parameter k. Secondly, RCP take also place during the homotopy process. 
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For stability reasons, some distances are removed and others are added on-the-fly according to 
numerical considerations. The elements of the reference of the CP are also modified during the 
homotopy process and in the following definition we put forward the reference as one of the 
characterizing elements of a RCP: 

Definition 3 (RCP) LetG = C[P,A] a PDSP. A RCP R of G is a quadruplet {I, C-, A+, Aq), 
where: 

• G- is a subset of G involving parameters A_, 

• A-|_ is a set of parameters called driving parameters, 

• /[P,A',Ao] is the CP o/C/[P,A'], 

. A' = A\A_UA+, 

• Ci\C are distance constraints called added constraints, 
and (I, C-, A_|_, Aq) is such that C_ = C \ Ci. 


In a RCP, the situation where a point is the intersection of two added constraints could 
not occur. This would create a new point which is not given in the initial statement. So, we 
will consider RCP that meet the following conditions: for each circle-circle intersection Pi+i = 
interCC(jpii ,SLi 2 ,pi 3 ,a.i 4 ), (i) 3^2 is not a driving parameter (i.e. 3 i 2 g A), and (ii) if a^i is a 
driving parameter {i.e. G A+), is a reference point {i.e. p^a G Aq) and is not used as a 
reference for another instruction. 

Given a RCP with a sole driving parameter in each instruction, it can always be modified to 
meet the conditions (i) and (ii) as follows: (i) is satisfied by rewriting instructions, and (ii) is 
satisfied by creating a new reference point for each driving parameter. A RCP for A 3 3 satisfying 
(i) and (ii) is R' = (/', {djst(p 3 , pe) = ag}, {k}, {pi, li, p(}) where /' is obtained from I by 
substituting P 3 = InterCC{pi,'k.,p 2 , ^ 2 ) by P 3 = /nterC'C'(p 2 , ag, p^, k). 

We focus here on the numerical step of the reparameterization method, that consists in finding 
values for A+ such that figures constructed by I fulfill constraints of C-. 

Let R = (J, C-, A_|_, Aq) be a RCP of G = G[P,A], A_ be the set of parameters of G_ 
and A' = (A \ A_) l±) A-|_, where l±) holds for a disjoint union. Given a branch b, we recall that 
[''l/(Ao,(A\ A_) W A+) is the evaluation of I on b. Both for the sake of readability and to 
make appear the different roles played by driving parameters and other parameters, we will note 
it ^^^/(Aq, A tt) A+) even if elements of A_ are not involved in I. When values Aq for Aq are 
explicitly fixed, we will note ^^h{A l±) A+) for ^^^/(Ao, A l±) A_|-). 

Let C- = {cj",..., c^}. We associate with the RGP R = (/, G_,A+,Ao) the numerical 
functions 

X K'" ^ 

/ cr(''''/(AWA+),A) 

A+,A ... 

V C-([''1J(AWA+),A) 



where 6 is a branch of / and numerical interpretations c~ are defined as in Eq. ([^. Since 
numerical interpretations c~ are G°°, functions ^^^R are C°° on the interiors of domains of 
definition of functions 

Given values A^° for A, there is a one to one cor respo ndence between real solutions of 


[b] 


R{A~^, A^°) = 0 for all branches b and real solutions of 
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3 Leading Homotopy by Reparameterization 

We aim at finding solutions oflying on the path S of|^to which belongs the sketch. Instead 
of using a path tracker to follow S in K™ x M, we propose to compute it indirectly by following 
a sequence of paths defined by homotopy functions constructed with a RCP. These paths are 
tracked in x M where d is the number of driving parameters of the RCP with d « m, what 
makes cheaper the path tracking. We give here a justification of this approach by showing that 
such paths are diffeomorphic to connected subsets of S. 

Sec. |3.1| enumerates assumptions that are required to make our approach valid. We define 
in Sec. |3.2| homotopy functions using RCP and characterize their domains of definition and 
boundary configurations in Sec. |3.3| We establish the link between their paths and S in Sec. |3.4[ 


3.1 Assumptions 

Let G = C[P, A] be a PDSP, and PI the homotopy function with interpolation function a. 
Let R = (/, C-, A+, Ao) be a RCP of G, where I = Ii is the instruction p 2 = 

interGL{pi, and for i > 2, Ii is the instruction Pi+i = interGC{pii,aLi2,pi3,aLn), where 

aj4 is either in A or in A^. Up to a re-indexing, the component ai of a interpolates values of 

Sj G A. 

The method presented here is valid under the following hypothesis on the interpolation func¬ 
tion a and the Jacobian matrix Jh of H. 


(hi) Jh has full rank on each point of 
(h2) supp+{a) is compact, 

(hS) Vi, j, ai{t) = aj{t) has a finite number of solutions, 

(h)) if t e supp+[a), ai(t) > 0, 

(h5) if i > 2 and Ii is s.t. G A, then ai 2 (t)aii(t) = 0 ^ ai 2 {t) ^ 0 or 0^4(t) ^ 0 for 
t G supp+{a), 

(h6) if i > 2 and R is s.t. G A+, then 0^2 (t) > 0 for t G supp+(a). 


In what follows, some instructions of the RCP will be changed in such a way that (hdj is satisfied 
only on a subset of supp+{a), and we will say that (h6) is satisfied on a given subset U of supp+{a) 
if (h6) holds for each t in U. 

Let us explain these hypothesis, (hi) and (h2) are the hypothesis of Cor.[^and they guarantee 
that paths ofj^are diffeomorphic to circles. When (hi) is satisfied, each point of II~^{0) admits 
a tangent and paths can be tracked with a numerical path-tracker. Here il(X, t) = F(X,a(t)) 
hence (hi) is strongly related to the rank of Jp, the Jacobian matrix of H. Since G is generically 
well-constrained, Jp has full rank on each point of F“^(0) for generic values of A (see Subsec. 2.1). 
Verifying that G is generically well-constrained and characterizing the interpolation function a to 
satisfy (hi) are both challenging problems and are beyond the scope of this paper. [12] justifies 
real homotopies thanks to the theorem of Sard. (h2) ensures that values t such that iJ(X, t) = 0 
has real solutions are in a compact interval (see Subsec. \2.2\. (h2) can be satisfied by setting 


a component Oi of a to —ct^ -\- -\- c)t -|- af^ with c > 0, and Uj for j ^ i to linear 

interpolations. 

Beside its influence on the topology of homotopy paths, a has an impact on the geometric 
configurations of the figures encountered in such paths. Here, we are using RCP to build nu¬ 
merical functions to track homotopy paths of H. The obtained functions are not defined in the 
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whole space of parameters, and the borders of their domains of definition are characterized in 
terms of geometric configurations, (hf), (h5) and (h6) restrain the geometrical configurations a 
path can pass trough. Configurations that can be encountered are detailed in Subsec. |3.3| (hS) 
is used to prove the termination of our algorithm. 

3.2 R-reduced Homotopy Functions 

Let R = (/,C-, A+, Aq) be a RCP and suppose values Aq for Aq are fixed. We call R-reduced 
homotopy the homotopies defined by 

['’liL«(A+,t) = ['’li?(A+,a(t)) = 0 {Hn) 

where 6 is a branch of R. The functions x M —>■ are called R-reduced homotopy 

functions and the connected components of solutions of IH a [ are called R-reduced paths. 

Considering Eq. ([^, the i-th component of is c“(^/(a(t) l±) A+), a(t)). Noting C 
X K the domain of definition of the function defined as 

A+,t [^'/(a(t) l±) A+) 

we state that the domain of definition of is and that is C°° in the interior of 

We are now interested in characterizing sets and their borders in terms of geometric 

configurations of objects constructed by I. A point (A+,t) £ x K belongs to the border of 
if its neighborhoods contain points of and points where is not defined. In general, 
is neither open nor close, and contains only a possibly empty subset of its border. We will 
call boundary of the subset of the border of that is in 

3.3 Domain of Definition and Boundary Configurations 

Since is the combination of functions we first characterize the domain of definition of 
the latter functions. 

Let p 2 = interCL{pi, ai, fi) be the instruction Ii. Since pi, li are part of the reference their 
values pith are fixed s.t. pi £ fi. Hence maps to the value ai of ai one of the intersections 
P 2 ,iA of li with a circle of radius Oi which center belongs to (see the configuration (cl) on 
fig. When ai > 0, there is an open neighborhood of oi where is C°°. Here we have 
oi = ai(t) and from assumption (hf), ai{t) > 0 on supp+{a). 

Let i > 2 and Pi+i = mterCC'(pii,aj 2 ,pj 3 ,ai 4 ) be the instruction R of I. maps to 

(pp, Oi2,pj3, 0 ^ 4 ) one of the two intersections p,p' of two circles. We make a disjunction on the 
number of intersections of the two circles. 

When the two circles are disjoint or coincident (with non zero radius), is not defined. 
When the two circles have exactly two intersections, is clearly C°° (see the configuration 
(c 2 ) on fig.|4l). 

Two configurations can lead the two circles to have exactly one intersection. The first one 
is when the latter circles are concentric with null radii, and 0^2 = 0^4 = 0. Recall that either 
an £ A or an £ A+. Suppose first an £ A. Hence an = an ft) and an = an ft), and assumption 
(h5) forbids the situation an ft) = an ft) = 0 when t £ suppn(a). Suppose now an £ A_(_. Hence 
an = an ft) and assumption (h6) forbids the situation an (t) = 0 while t is in a subset for which 
it holds. As a consequence, an = an = 0 does not happen when (h5) and (h6) hold. 

The second configuration is when the two circles are tangent, and at least one circle has a 
strictly positive radius. The two centers of circles and their intersection are collinear (see the 
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Figure 4: Three geometric configurations for InterCL and InterCC instructions. 


configuration (c3) on fig. . Clearly, this configuration characterizes the boundary of the domain 
of definition of the mapping , aj 2 , p^s, aj4) and is called a boundary configuration of fi. 

Since I is the combination of its instructions, a point {A+,t) is in the boundary of only 
if a boundary configuration holds for (24_|_,a(t)) for at least one instruction li with i > 2. We 
will say in this case that the figure ^^^I{a{t) l±) ^4+) presents a boundary configuration of I or R, 
or that {Aj^,t) leads to a boundary configuration of / or R. 


3.4 R-Reduced Paths 

The point here is to characterize i?-reduced paths, and to link them with paths of|^ To achieve 
this, let us define the mappings x K —>• K™ x K as 




and (f' : M"* x K —>• x K as 




/ c^(X) \ 

4(X) 

V t / 


(3) 


(4) 


where c^,...,cj| are the d added constraints of R. are C°° on and p' is C°° on 

X M. Consider the following remark, that is a consequence of the characterization of a 
boundary configuration. 


Remark 1 Let (X,t) be a point of P[ ^(0) and R a RCP. X does not present any boundary 
configuration if and only if it exists a unique branch b and a unique point (A_|_,t) = p'{X,t) s.t. 
{X,t) = ^^^(4+,^). X presents a boundary configuration if and only if there exist at least two 
branches bl and b2 and a unique point {A^,t) = p'{X, t) s.t. {X, t) = t) = t). 

Let us justify Rem. If {X, t) is a point of (0), (^4+, t) = p'{X, t) is unique by definition 
of p'. Suppose X does not present any boundary configurations: each point of X is constructed 
by R by intersecting two circles having two different intersections. If 51 52 are two branches 

of R, they correspond to different choices of intersections and ^ '^^'^^p{A+,t). If X 

has a boundary configuration, at least one point p of X is the intersection of two circles in 
configuration (c3). The two branches 51 52 corresponding to the same choices for each point 

but for p are such that {X,t) = ''^^^p{A^,t) = t). 

We show now that i?-reduced paths are locally diffeomorphic to paths of[^ then we extend 
the latter diffeomorphism to a global diffeomorphism between pieces of i?-reduced paths and 
pieces of paths of|^ 
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Lemma 1 Let {A^,t) be s.t. ^^^Hji{A^,t) = 0, and X = l±) A+). If{A~^,t) does not lead 

to a boundary configuration of R, it exists a neighborhood U of (A"*", t) and a neighborhood V of 
{X,t) such that CiU is diffeomorphic to H~^{0) n V. 

Proof of Lem. 0 Let X = l±) A+). Then H^{A^, t) = 0 => Pl{X,t) = 0, and it exists 

a homotopy path S of|^to which belongs {X,t). Let "H- be the system of equations having all 
the equations of but the ones corresponding to constraints of C-, and S- be the set of its 
solutions. S- is a d + 1-dimensional smooth manifold, and {X, f) € S- and S C S- hold. As a 
consequence, 5 is a 1-dimensional smooth submanifold of S-, and in any open neighborhood of 
{X,t) in iS_, 5 is a 1-dimensional smooth manifold. 

Since (A+, t) does not lead to a boundary conhguration, is C°° on an open neighborhood 
C X M of (A+,t). Let V = Clearly V C and (A, t) S V hold. We show that V 

is an open neighborhood of (A, t) in S-. Let —?► x K be the restriction of p' to S-. 

tp' is C°° on S- and p'iV) =14. Hence V is an open neighborhood of (A, t) in as the inverse 
image of an open neighborhood. 

We finish the proof by remarking that the mapping ■ 14 V, which inverse is , is 

a diffeomorphism. □ 

The following Prop, extends the property of Lem. [^to a global property. 

Proposition 1 Let S be a homotopy path oM and R be a RCP. Let S' S be a connected subset 
of S that does not contain any figure with a boundary configuration of R. Then it exists a unique 
branch b, a unique R-reduced path and a subset C ^’’^S such that is diffeomorphic 
to S' by 


Proof of Prop. Let (A, t) be a point of S' . Since A does not present a boundary configuration 
of R, it exists (see Rem. a unique branch b and a unique point (A+,t) = p'{X,t) such that 
= (A, t). Hence H R{A.^. , t) = 0 and (A_|_,t) belongs to an homotopy path of 

Since S' is connected and p' is C°°, p'{S') is a connected subset of x M. Moreover, p'{S') 
belongs to the interior of otherwise a point of p'{S') would belong to the boundary of 

and would lead to a boundary configuration of R. As a consequence, '^^^p is well defined on 


p'{S'). 

We show now that p'{S') C If it is not the case, it exists at least another branch 62 
and another i?-reduced path with p'{S') C ([^^5 U and n ^ 0 holds from 

the local diffeomorphism property. Let {A'_^,t') G (^^^5 n Then {A'_^,t') is such that 

{X',t') = ^^^p{A'^,t') = ^^^^p{A'_^_,t') and {X',t') G S' presents a boundary configuration of R 
from Rem. [T] 

Let us defined ^^^S' as p'{S'). It is a connected subset of ^^^5, and = S'. Finally 

^^^p is injective and is a global diffeomorphism from ^^^S' to S'. □ 

Prop. [2 states that it is possible to compute a part of S that does not contain any boundary 
configurations by following a path of 

Now, assuming that the figures of S presenting a boundary configuration of R are in a finite 
number (boundary configurations of type (c3) can be each described by a polynomial equation, 
hence figures presenting a boundary configuration of R can be seen as the solutions of systems 
of m -I- 1 polynomials involving m -|- 1 variables) the set S can be written 5^ U (A^,t^) U 5^ U 
... U 5" U (A", t"), where (A®, P) are such that A® presents a boundary configuration of I, and 
S^ Q S are connected, pairwise disjoint and does not contain any figure presenting a boundary 
configuration. From Prop. [^each path 5* can be computed by following a path ^^*^5 of 
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Figure 5: Left: the sequence of i?-reduced paths corresponding to a path S of|^for Black 
circles mark points leading to a a boundary configuration. Right: geometric configurations 
corresponding to points (a), (b), (c). 


hence S can be computed by following the sequence of paths ^^*^5 and making the appropriate 
branch changing in points 

The left part of fig. shows the sequence of i?-reduced paths corresponding to a path 
of 1^ for the PDSP -K 3 3 . Its right part shows geometric configurations near a point leading to a 
boundary configuration. 


4 On-The-Fly Change of the RCP 


The method roughly depicted above hides a pitfall: it leads to follow i?-reduced paths until 
a boundary configuration is reached. But i?-reduced paths are numerically bad behaving near 
boundary configurations. 

It can be easily seen by considering a function that maps to a positive real number k the 
positive ^-coordinate of the intersections of the unit circle with a circle centered in (0, 2) of 

lOfc^ — 9 

radius equal to k. Noting y this function, we have y{k) = -- and y'(k) = 

— 4^3 _j_ 2Qk 

The domain of definition of y is [1,3] which boundaries 1,3 corresponds to 


sy^FTToF - 9 

configurations where the circles are tangent; we have lim]~^^i 3^y'(k) = + 00 , and it is not due 
to the chosen system of coordinates. 

Such unbounded values of derivatives highly affect the efficiency of a numerical path tracking, 
that proceeds by approximating a path by its tangent. 

We propose here to introduce a measure of the distance from a figure of 5, or from a point of 
a i?-reduced path, to a boundary configuration, and to stop the tracking process of a i?-reduced 
path when this distance is smaller than a real parameter a. Then we change either the RCP or 
the values of its references in a way that the distance to a boundary configuration is greater than 
a. The new R-reduced homotopy path is then followed to compute the path S. This process is 
repeated each time the distance to a boundary configuration is smaller than a. 

The distance to a boundary configuration is defined in Sec. |4.2[ and our algorithm to change 
on-the-fly the RCP is described in Sec. |4.3| We first give an intuition of our approach on the 


example K3 3 in Sec. 4.1 
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Figure 6: Center: a zoomed view of a sequence of i?-reduced paths for K 3 3 . Pieces in solid lines 
are followed by our algorithm. Black circles mark points leading to a a boundary configuration. 
Left and Right: geometric configurations corresponding to points {a),...{d). Dashed double 
arrows correspond to constraints removed by our algorithm. 


4.1 Overview of our Method on an Example 

The middle part of Fig. shows a zoomed view of pieces of i?-reduced paths {R is given in 
Sec. |2.3.3 | for 1^3,3. When following this union of paths from the sketch, the point (a) is first 
reached. The points P5,p4,p2 in the figure constructed in (a) is shown in the right part of 
Fig. § The two circles constructed when evaluating the instruction ps = interCC(p 2 , ag, p4, a4) 
are almost tangent, and is “too close” to a boundary configuration (this notion will be detailed 
in Sec. |4.2[ ). 

To avoid the point (b) that leads to a boundary configuration, a new driving parameter and 
a new reference point is added to the RCP: the instruction ps = mterC'C(p2, ag, p4, a4) is re¬ 
placed by Ps = mterCC'(p2, ag, P4, k2), where P4 is a new reference point and k2 is a new 
driving parameter. The constraint distance(p4, ps) = a4 is added to the set of removed con¬ 
straints to guarantee that it is fulfilled by constructed figures. For an appropriated placement 


of P4, presented in Sec. 4.3.3 the figure constructed by the CP is not too close to a bound¬ 
ary configuration. The new RCP defines new paths, that can be followed while avoiding the 
boundary configuration of point (6) (see right part of Fig. |^. When reaching the point (c), the 
original instruction ps = jnter(7C(p2, ag, p4, a4) can be restored while staying “far away” from 
a boundary configuration. The piece of path between (a) and (c) in the middle part of Fig. 
is drawn in dashed line to underline that it is a projection of the path that is followed with our 
algorithm. 

Suppose now that the point (d) (see central part of Fig.[^ is reached. The constructed figure 
is too close to a boundary configuration of the instruction pg = interCC{p 2 ,a 2 , Pi, where 
k is a driving parameter. In that case, the point p^ is moved in order to avoid the boundary 
configuration of point (e) (see left part of Fig. [^. The piece of path after (d) is drawn in dashed 
line to figure out that it is no longer the path that is followed by our algorithm. 


4.2 Distance to a Boundary Configuration 

Let li be Pi+i = interCC{pii,ai 2 ,pi 3 ,aLii). We associate with R the real function 7^ taking its 
values in [0,1]: 

dist{pi+i,{pii,pi3)) 

lii.X) “ - -- ■■■ , 

maa:(pi+ipp, p*+iPi3) 
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Algorithm 1 Path tracking with RCP swapping 

Input: A RCP R = (I, C-, A-|_, Aq), an interpolation function a, a sketch a. 
Output: a list of solutions Lsoi 
1: Let t = 0, A = A®'' 

2: Find A+, Ao, h s.t. a{t) l±) A+) = A, fix Ag to Ag 

3: while True do 

4: Follow the path C from (A+,t) while checking 

5: if passes trough the sketch then 

6 : return Csoi 

7: else if passes trough hyperplane t = 1 then 

8: Append current figure to Cgoi 

9: Set (A+,t) to current point 

10: else if current point (A™’',t™’') € satisfies (scl) or ^*’^7/(A™’", < a then 

11 : Let A = tt) A™’') 

12: Apply Algo. to change R = (/, C-, A+, Ag), and obtain A+, Ag 

13: Find b s.t. ^^T(Ag, a{t) l±) A+) = A, fix Ag to Ag 


where dist(jpi+i, (ppjPis)) is the distance between Pi+i and the line (pp,Pi 3 ). 

7 i is defined and continuous on A when points {pi,Pi^,Pi3} C A are not coincident. Notice 
it never happens when R is an instruction of a RCP for which assumptions (h5) and (h6) hold. 
7 i vanishes only on figures presenting a boundary configuration of R. 

We associate to the function 7 / that measures the distance of a hgure A to a boundary 

configuration of / defined as 

77(X) := Tnin2<i<ni(X). 

7 / is defined and continuous at least when assumptions (h5) and (h6) hold, and vanishes on 
figures presenting a boundary configuration of I. In the following, we will note ^^^ 7 (A+,t) for 
7 /(''’'/(a(t)WA+)). 

4.3 Path Tracking with On-The-Fly Change of the RCP. 

Algo, [^describes the main process of our method. We consider a PDSP G = (^[X, A], a RCP 
R, a sketch A®^ and an interpolation function a from A®^ to A®°. We assume that assumptions 
(hi) to (h6) are satisfied for a and R. 

In the step 2, value Ag for elements of Ag are read on A and fixed. Values A_|_ for A+ are 
obtained by evaluating ip'{X,t). In steps 2 and 13, the branch b is found by evaluating one by 
one instructions R on each branch, and keeping for each bi the choice that leads to construct A. 

When entering in the step 4, a point (A, t) of the path S of|^to which belongs the sketch is 
known as well as a point (A+,t) with = (A, t). We temporary assume (h7): X is not 

closer than a to a boundary configuration. 

In the step 4, an abstract path-tracker is used to follow in a given orientation the i?-reduced 
path s.t. (A_|_,t) £ We assume that it allows to compute in a finite number of itera¬ 
tions the connected subset ^^^5' C [^1(5 s.t. points of ^^^5' are not closer than a to a boundary 
configuration, and that it stops when a point at a distance a to a boundary configuration is 
reached or when the stopping condition (scl) described below is satisfied. It is also assumed that 
it is possible to detect when passes trough hyperplanes t = 0 and t = 1 , and to get exact 
intersections with latter hyperplanes. 
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Algorithm 2 Change RCP 


Input: A RCP (/, C-, A-|_, Aq), a table T of original instructions, a current figure X, a 
Output: New values A+, Aq 
1 : for 2 < i < I do 


2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 
23 


Let li G I he Pi+i = mterC'C'(pii,aj2,pj3,aj4) 
if 7 i(A) < a and T[i] == 0 then 
if SLii G A-f. then 

Apply Algo. a to obtain 0^4, 
else 

Let p'3 be point and a^+i a length parameter 
Let p'3 = Pi3 and aa+i = 0^4 
Let /' be Pi+i = mterC'C'(pii,aj2,p'3,ad+i) 
Apply Algo. 3 with f and /' as inputs 
Apply Algo. 4 to obtain ad+i,p'^3 
else if not T[i] == 0 then 

Let /' = T[i] associated with 7 ' 

Let /' be p^+i = mterCC'(pii,a,;2,p'3,a'4) 
if 7i(A) > a then 

Apply Algo. with f and /' as inputs 
else if 0^2 < a'4 then 

Apply Algo. i with li and /' as inputs 
goto step 3 

else if 7i(A) < o then 

Apply Algo. a to obtain 0^4,pj 3 
Actualize A^,Aq and d = |A+| 

return A+, Aq 


[> assume a ,2 > a ,4 


\> restore /' 


> change driving parameter 


t> move point 


The tracking process stops when a point (A+, f) of is at a distance a to a boundary 
configuration. In this case a new RCP, or at least new values Aq for reference points and 
driving parameters, is computed thanks to Algo. described in Sec |4.3.1[ and the process re¬ 
enters in step 4 with a new RCP R, new values Aq,A^ and a new branch b such that (A_|_,t) 
is not closer than a to a boundary configuration. Hence assumption (h7) is satisfied when re¬ 
entering step 4. If (hi) is not satisfied when performing for the first time step 4, Algo. is 
directly applied. 


Stopping condition (scl) is detailed in Sec. |4.3.2 It is satisfied when it is not possible to 
ensure that assumption (h6) holds. When (scl) is satisfied, Algo. changes the RCP in such a 
way (h6) holds on the computed path. 


Sought solutions are found when passes through the hyperplane t = 1, and the overall 
process is stopped when a point (A_|_,t) of is s.t. ^^^I{a{t) W A+) = Latter termination 
criterion is checked each time the the hyperplane t = 0 is crossed. 


The orientation used to follow the new i?-reduced path is chosen in order to avoid backtrack on 
S. We will discuss practical details of path-tracking in Sec. 7^ We now focus on the description 
of the way the RCP is changed. 
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Algorithm 3 Swap instructions 

Input: A RCP R = (/,C-, A+, Aq), an index i, a new instruction a table T of original 
instructions 

1 : Let li G I he Pi+i = interCC{pii, Sup, p^s, 8144 .) 

2 : Let /' be Pi+i = mterC'C(pii,aj2,p'3,a'4) 

3: Let c be the constraint distancelpi^iypp) = aj4 

4: if not c G C- then > introduce new instruction 

5: A+^A+U{a'4} 

6 : Ao G- Aq U {p^a} 

7: / ^ (/i U .. . U h-i U /' U /*+i U . . . U h) 

8: C. ^ C_ U {c} 

9 : T[i] G- li 

10 : else > restore original instruction 

11 : A+ A+ \ {a' 4 } 

12 : Ao^AoMp's} 

13: / ^ (/i U ... U /,_i U /' U /,+i U .. . U Ii) 

14: {c} 

15: T[i] G- 0 

16: return 


4 . 3.1 Changing RCP 

In Algo[^ when a point {A+,t) s.t. = a is reached, the RCP or at least the values 

of the reference are changed. The basic principle of the mechanism that changes the RCP or the 
reference is to identify the instruction(s) R s.t. ji(X) = a, where X = l±) A+). 

Let R be Pi+i = interCC{pii, 8 .^ 2 , p^ 3 ,a.i^) s.t. 7 i(A) = a. Then either aj 4 G A_(. and 
Pi3 G Aq is a reference point involved only in li, or aj4 G A, Pj3 ^ Aq and we suppose without 
loss of generality that 0^2 > 0^4 (otherwise arguments of the instruction are swapped). 

In the first case, new values for pjs, aj 4 are computed s.t. 7 i(A) > a thanks to Algo.|^ In the 
second case, R is exchanged with the instruction /' defined as Pi+i = interCC{pii , 812 , p' 3 , a^+i) 
and involving the new driving parameter ad+i and the new reference point p'a. Then values 
ad+i^p'i 3 are computed s.t. 7 i(A) > a, where 7 ^ is the distance to a boundary configuration 
of This relaxation is counterbalanced by adding distance{pi+i,Pi 3 ) = aj 4 to the set C_ of 
removed constraints of the new RCP. 

A table T, that does not appear in Algo. to ease its description, is used to save the 
instructions of the original RCP {i.e. given as input of Algo. [^. Entries of T are initially 
empty, and each time an instruction R of the original RCP is exchanged with /(, R is stored in 
the i-th entry of T. When R could be restored while ensuring that the distance to a boundary 
configuration stays greater than a, li is restored and T\i] is re-set to 0. Hence an instruction R 
of the current RCP is an instruction of the original RCP if T\i] is empty. 

Algo.details the mechanism to change a RCP, and Algo, [^details the way instructions are 
swapped. They both modify in place the RCP. Notice that it could exist several indices i such 
that 7 i(A') = a, and that Algo. is designed to take it into account. 


4 . 3.2 Stopping Condition (scl) 

The tracking process in Algo.[^also stops when the condition (scl) is satisfied for a current point 
(A“’',t™’'). We define here this stopping condition. 
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Algorithm 4 Shift Reference 

Input: An instruction Pi+i = mterC'C'(pii , ai 2 , p^s, aj 4 ), values Pi+i,pii,ai 2 ,pi 3 ,ai 4 for 
Pi+l) pp ) , aj4 . 

Output: New values ^^ 3 , 0^4 for 
1 : Let p be the projection of Pi+i on the line (piijPis) 

2 : Let p '3 = Pi+i + ai 2 
3 : Let a'4 = 0^2 
4: return p' 3 , 0^4 


Pi+i 



Figure 7: New placement of points given by Algo. p (resp. p') is the projection of Pi+i on the 
line (piijpp) (resp. (ppjp'a)) and v (resp. v') is the distance from to p (resp. p'). 


Let li be an instruction of I s.t. T[i] 0. Hence f has been introduced by Algo. to 
replace the original instruction I[ = T[i]. Let li be Pi+i = mterC'C(pp,aj 2 ,p' 3 ,a' 4 ), and /' be 
Pi+i = *nterC'C(pp,ai 2 ,pj 3 ,aj 4 ). Recall that in our homotopy context, values for ai 2 ,aj 4 are 
Oitif) and 0 ^ 2 ( 1 ). 

(scl) is satisfied if it exists 2 < i < I s.t. T[i] f 0 and 0 ^ 4 (t“’’) > 0 ^ 2 ( 1 ™’’). When 
(scl) is satisfied on a point of a path, Algo. is called. Unless the instruction li making 

(scl) to be satisfied has been restored in steps 15-16, the if condition in step 17 of Algo. 
is satisfied, and steps 18-19 are performed: the original instruction /' is restored, and when 
entering step 6 , its arguments are swapped {i.e. Pi+i = interCC{pii,ai 2 ,pi 3 , 8 ^ 4 ) is replaced by 
Pi+i = interCC{pi 3 ,ai 4 ,Pii,ai 2 )). Then a new driving parameter is introduced. 

When returning to Algo. [Rafter Algo, [^have been performed, the stopping condition (scl) 
is not satisfied. We will state in Sec. |5.1|that this mechanism ensures that (h 6 ) holds. 


4.3.3 Shifting Reference 

Algo. a is called in steps 5, 11 and 21 of Algo.[^ It computes values for driving parameters and 
reference point of an instruction R in order that the constructed figure X satisfies 7 i(A) > a. 
The following proposition states that this goal is achieved after applying Algo. if 0 < a < 

Proposition 2 Let 7 ^ be the distance to a boundary configuration associated with the instruction 
Pi+i = interCC{pp, ai 2 ,pi 3 ,ai 4 ) and (pi+i,pii, ai2,pi3, 0 * 4 ) be values s.t. jt{pt+i,pp Pis) < a. 
If 0 < a < ^ and p '3 has been obtained with Algo. 0 then 7i(pj+i,pp,p'3) > a. 

The proof of Prop, [^is depicted in fig. where p is the projection of Pi+i on the line (pp,) 
and V is the distance from Pi+i to p. Let p '3 be the new placement of Pi 3 obtained with Algo. Hi 
The criterion 7 ^ is computed by considering the distance v' from p^+i to p' where p' is the 
projection of pi+i on the line (pp,p' 3 ). One has i + 2 aT result follows. □ 
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5 Correctness and Termination 


Let S be the homotopy path of to which belongs We show here that Algo. 

terminates, and that all solutions {X,t) of with t = 1 lying on S are in Cgoi at the end of 
Algo. 

Algo. involves an iterative path tracker that is assumed to track in a finite number of 
steps a i?-reduced path if it is a manifold, and if its points are not closer than a to a boundary 
configuration. The correctness and the termination of our method is proved when considering 
such an abstract path tracker. Sec. [^proposes an implementation of such a path-tracker. 

Let R = (/, C-, A+, Aq) with I = {Ii)i<i<i be the RCP given as input of Algo. The 
latter procedure computes a sequence {R^)j^N of RCP where A C N*, = i? and R^ = 

(A, C:L, A-^, Aq) with P = A sequence of connected pieces of R-^-reduced 

paths is followed, and we note bj the branch of A such that = 0. We will note 7/j the 

distance to a boundary configuration associated with R and for q/j (a(t) tt) 

A-^)). We will note and (/?' the mappings defined in Eqs. and [^specialized to the RCP 
R-'. Let finally TTt be the projection with respect to the t-coordinate. 

The main points of the proof are: 

(z) is well defined on , 

(ii) if G and a < | then (A-^,t) > a, 

(in) {S^)j^N are 1-dimensional manifolds, 

(v) A is a finite subset of N and UjeAf 


Remark that the sto pping condition (scl) is not satisfied on a point {A\,t) G , for j G A. We 
will prove (z) in Sec. 5.1 by showing that (h^) and (h5) hold for R-^, and (h6) holds for R^ on 
7rt(5'’). Then (zz) holds thanks to Prop.[^ 

(zzz) and (zu) are consequences of Prop.[^ is a diffeomorphism from to a connected 

subset of S. Notice that (zz) and (zzz) are the two conditions under which the abstract path- 
tracker used in Algo, [^computes SP 

{v) is proved in Sec. |5.2[ It has as a direct consequence that Algo. terminates and all 
solutions (A, <) ofwith t = 1 lying on S are found. 


5.1 Proof of Point (i) 

As stated in Sec. |4.2[ the distance to a boundary configuration associated with a RCP is well 
defined when (h4), (h5) and (h6) hold. It is established in the following proposition, and point 
(zz) follows as a corollary. 

Proposition 3 If assumptions (h4), (h5) and (h6) hold for R, thenMj G A, (h4) and (h5) hold 
for R^, and (h6) holds for R^ at least on 7rt(5^). 

Proof of Prop. Let j G A. The first instruction of the RCP is never changed in Algo. 

= Ii follows and assumptions (h4) holds for R^. 

Let 2 <i < I and If be the instruction Pz+i = znterC'C'(pp, aj 2 , p^a, aj4 ). If aj4 ^ A^ then 
If = Ii- Since (h5) holds for R, (h5) holds for Rj. 
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Suppose now 3^4 G If acii G A+, i.e. aj 4 is a driving parameter of R, then If = li and 
0^2 (t) > 0 since (h6) holds for R when t G supp+{a). 

Otherwise, If is not an original instruction. Let Pi+i = interCC{pii, ai 2 ,p' 3 ,a- 4 ) be the 
original instruction R and 0^4 be the interpolation function for values of a' 4 . 0^2 (t) and 0^4 (t) 
does not both vanish according to (h5). Since condition (scl) is not satisfied for {A^^,t) G SR 
0^2 (t) > 0^4 (t) holds and 0^2 (t) does not vanish on 7 rt( 5 '^ ). Thus (h6) holds for Rj on Trt{S^). □ 

Corollary 2 If assumptions (hf), (h5) and (h6) hold for R, thenMj G N, is well defined 

on . 

5.2 Proof of Point (u) 

We consider first the case where A = {1,..., n}. The termination condition of step 6 of Algo. 
is reached, hence the set UjeAf diffeomorphic to a circle. Since IJ^-g ^ C 5, 

UjeAT = S follows. 

We consider now the case A = N* and we show that it never happens. Algo. constructs 
a sequence {{AR,P))j^]\f of points s.t. either {AR,P) = a or (scl) is satished, and 

,P~^) > a and (scl) is not satisfied. Consider the sequence {{XRtR)j^N where 
{XRtR = From point {iv), {{XRP))j^f^ is a sequence of points of S. 

From Cor. 5 is diffeomorphic to a circle, hence 5\ {(A'*^,0)} is diffeomorphic to a bounded 
open interval, and it exists a diffeomorphism S :]0,1[—>■ S \ {(A'^^,0)} that maps to s e]0,1[ a 
point of 5 \ {(A'’^,0)}. Reciprocally, S~^ maps to a point (XRP) a real number s^ g]0, 1[. 

We show that the sequence {sRj^n satisfies Vj > 2, and does not have any accu¬ 

mulation point. As a consequence, it cannot be inhnite. 

Suppose it exists j > 2 s.t. s^ = s^~^, hence (XRp) = {X^~^,P~^). From Cor. 7/3 is 

well defined on X^ = X^~^ and 77 j(A-^“^) > a hence (scl) is not satisfied for (A-^”^, and 

either 7 / 3 (A-^) = a or (scl) is satisfied for (XRp) hence a contradiction follows. In Algo. 
paths S^ are followed with an orientation that ensures a progression along S, hence we have 
s^ > sAi, 

Suppose now that the sequence {sRj^n has an accumulation point s*, hence {{XRp))j^fq 
has an accumulation point (A*,t*), and it exists a subsequence {{XRPRj^n* , with A* C A, 
converging to {X*,t*). From assumption (hS), there is an index s.t. Vj > j^,Vl < < m 

the sign of aii{P) — aii{P) does not change. Hence for j > j^, for each instruction = 
inter CCfpp^ 3^2, of P with i>2, aR > aR and since assumptions (h5) and (h6) hold, 

it exists r > 0 s.t. Vj > j^,Vi > 2,af2 > r. Now, for each e > 0, it exists an index j^ > j^ s.t. 
Vj > j^, points and distances between points vary no more than e between A^ and X^~^. Remark 
that since signs of aii{P) — ai 2 {P) does not change when j > j^ grows, (scl) is satished neither 
for A-l nor for A^i and it follows that 7/3 (A-l”^) > ^ (from Prop. and 7/3 (A-1) = a < ^ 
when j > 2. Taking e sufficiently small (strictly less than (^ — Q;)r) leMS to a contradiction. □ 


6 Generalization to 3D PDSP 

3D PDSP Ht well to the method depicted in this paper: results of [3] as well as reparameterization 
approach stay valid. Given a PDSP G in a 3D geometric universe, solutions of G up to rigid 
motions are found by hxing a reference consisting in a point pi, a line R and a plane pR. Values 
pi,li,pli are hxed s.t. pi G R and h G pR. A CP of G has the structure: 

(R) P2 = inter S'L(pi,ai,R) 
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Figure 8 : Placement of point p '5 when swapping a interCS instruction. Dashed lines does not 
belong to the plane of circle Ci+i. 

{h) P 3 = interS'S'P(pi, a 2 , P 2 , as, pl^) 

{h) P4 = znterS'5'5'(pi,a4,p2,a5,p3,a6) 


Pz+1 ^nter 5**5^iS^(pp, aj2 , Pj3 , a^4 , p^s, a^6 ) 


where interSL is a sphere-line intersection, interSSP is the intersection of two spheres and one 
plane, and interSSS is a three spheres intersection. 

If li is an interSSS instructions, it is decomposed into the two instructions 


(Ili) Ci+i = interSS{pii,a.,2,p,3,ai4,) 


{I2i) pi+i = interCS{c{ci+i),r{c,+i),Pi5,ai6) 

where Ci+i is a circle, interSS is a sphere-sphere intersection, c(ci+i) and r(ci_|_i) are respectively 
the center and the radius of c^+i, and interCS is a sphere-circle intersection. 

interSSP and interSS can be seen as a 3D extension of an interCC instruction in the 
2D case. Hence when assumptions (hi) to (h6) hold, boundary configurations encountered on a 
homotopy path S of|^for these instructions are the same than boundary configurations in the 2D 
case. The function that measures the distance to a boundary configuration of such instructions is 
the natural extension of the one associated with an interCC instruction. When such instructions 
are swapped, reference points are fixed as in the 2D case. 

Consider now an interCS instruction 72^. A boundary configuration is reached when the cir¬ 
cle and the sphere are tangent. The radius of the center is null only for boundary configurations 

of Hi- We associate with (72^) the function 7^ defined as 7i(X) = --- P d il? --- 

max{p^+ic(a,+i), Pz+iPi5) 

where p' is the projection of p^s on the plane to which belongs c^+i, and p is the projection of Pi+i 
on the line passing by c(cj+i) and p'. 7 ^ is defined on figures that are not a boundary configura¬ 
tion of 71i. If a^B is not already a driving parameter, Pi+i = mterCS'(c(ci+i), r(cis_i), p^s, a^e) 
can be swapped with Pi+i = interC'5'(c(ci+i),r(ci_|_i),p'g,a'6), and the value p'g for the new 
reference point p'g is set as p'g = where m is the greatest value between r(ci_|_i) and 

QiB, and a'e is set as m, as illustrated in Fig^ It is then easy to state a proposition equivalent 
to Prop. for this placement of point. 
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7 Implementation and results 

Our method has been implemented in C++, giving rise to a program that accepts a PDSP G, 
a RCP, a sketch and provides solutions of G. The path-tracking is achieved by a prediction- 
correction method with an adaptive prediction step that is described in Sec. |7.1[ It requires 
to compute partial derivatives of the function which appears to be one of the most time 

consuming step of our method. We propose in Sec. |7.2| to exploit the acyclic nature of a CP to 
optimize this operation. Numerical results related to the solving of four PDSP in 2D and 3D are 
given in Subsec. |7.3| It confirms the efficiency of the approach of [9] to provide several solutions 
(sometimes all of them) of problems that resist to divide and conquer methods and are too large 
to be solved by classical numeric solvers providing all the solutions. Using a RCP as proposed 
here brings an important speed-up of this approach. 


7.1 Path tracking 

Homotopy and i?-reduced paths are followed thanks to a classical prediction-correction method: 
prediction is performed along the tangent of the path by an Euler predictor with step 6 S 
[5min,dmax] and correction by Newton-Raphson iterations. The step 6 is doubled (respectively 
halved) if 2 S < Smax (resp. 6 > 2 Smin) and if the previous correction step did succeed (resp. fail). 
The Jacobian matrices that are required both in prediction and correction steps are numerically 
computed with finite differences. 

In Algo. when entering for the first time in the main while loop, an orientation (ie. one 
of the two unit vectors of the tangent) to follow the first i?-reduced path is arbitrarily chosen. 
When entering in the while loop after the RCP has been changed for the j -I- 1-th time, the 
orientation has to ensure the progression along S. To determine the appropriated orientation, 
the last unit vector of the tangent used to track is “translated” in the new space where 5'^+^ 
is tracked with the application o with notations of Sec. 

Notice that this simple path tracking algorithm does not avoid jumps between paths. Ap¬ 
proaches using interval arithmetic (see [101 [3 [T3]) could be used to certify the path tracking. 
Here we suppose that Smax is small enough to follow considered curves while avoiding such jumps. 


7.2 Differentiation of the CP 

When tracking a i?-reduced path, most of the computation time is spent in the evaluation of the 
underlying RCP. Most evaluations intervene in the computation of Jacobian matrices by finite 
differences that needs about d evaluations, where d is the number of driving parameters. Such 
matrices are computed at each prediction step and at each iteration of the Newton-Raphson 
method in a correction step. Here we exploit the acyclic computation scheme of a RCP to 
improve its evaluation. 

Suppose driving parameters = (ki,..., k^;) of R appear in instructions Ip,..., I^d with 
i‘^ >...> A. When computing with finite differences the derivative with respect to kj of the 
numerical function associated to R, the geometric objects of the figures resulting of the two 
evaluations differ only if they are produced by instructions R with i > P since lid is the first 
step involving k_, . 

Hence a manner of optimizing the differentiation of a RCP is to evaluate it entirely a first 
time and then to compute the partial derivatives with respect to kj by evaluating the RCP from 
the step U, for j from d to 1. Our implementation incorporates this optimization. 
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Figure 9: Octahedron (left) and disulfide (right) problems. Edges are distance constraints. 

7.3 Results 

We give here numerical results concerning the solving of four PDSP, one in 2D and three in 
3D. The method depicted here consists in computing the path S of to which belongs the 
sketch by using a RCP to track S in the space of driving parameters instead of tracking it in 
the space of all coordinates. These two approaches (with and without RCP) yield the same 
number of solutions. However, using a RCP brings an important gain in term of running times 
as it appears in our experiments. For each problem we also give the running time and the 
number of solutions obtained when solving the system with a classical homotopy method 
implemented by the free software HOM4PS-2.0 (see [H])- We did chose homotopy solving as a 
witness method because as far as we know, it is the sole approach allowing to find all the solutions 
of large undecomposable problems. We did choose HOM4PS-2.0 to implement it because among 
other free softwares implementing homotopy, it seems to be faster to solve sparse systems of 
polynomials. 

Be given a real solution of a PDSP, the elements of its orbit by the action of the group of 
reflections through x and y axis (resp. (x,y), {y,z) and (x,z) planes in the 3D case) are also 
solutions of the PDSP. In our experiments, the solutions found on S belong to different orbits. 


7.3.1 Problems and parameters settings 


The goal of the octahedron problem is to construct a solid with 6 vertices, 12 edges and 8 
triangular faces knowing the lengths of its 12 edges (see Fig. [^. This problem is used in [3] and 
is related to the parallel robot called Gough-Stewart platform. It results in a system p^sol of 12 
equations with 12 unknowns. 

The second problem comes from molecular chemistry and is picked up from m- Coordinates 
of 8 points in the 3D space have to be found knowing 18 distances. It corresponds to a disulfide 
molecule (see right part of Fig. [^. A valuation of parameters is exhibited in m that leads to 
18 solutions up to reflexions all found by a bisection method in more than 10 minutes in HZl. 

Dodecagon and Icosahedron problems are illustrated on Fig. 
system Ea with 21 equations. The system EJ 


10 The former gives rise to a 


associated to the icosahedron problem involves 

30 equations. 

Values of the sketches and of parameters are given in appendix [a| 

Interpolation functions have been chosen such that am{t) = + (a^ — + 2)t + and 

ai(t) = (1 — t)al^ + taf°, for 1 < i < to. The value for a has been set for each problem to 0.1, 
what seems to fit well to our algorithm. The prediction steps vary in the interval [1~^°, (5maa;]- 
As stated above, Smax has to be chosen small enough to avoid jumping between different paths. 
The values Smax = 0.1 when tracking S without RCP and Smax = 0.05 when using a RCP have 
been chosen after several trials. Further details and discussions concerning the prediction steps 
are given in|7.3.3[ 
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Figure 10: Dodecagon (left) and icosahedron (right) problems. Edges are distance constraints. 
Table 1: Sequential running times on an Intel(R) Core(TM) i7-5600U CPU @ 2.60GHz. 



Octahedron 

Disulfide 

Dodecagon 

Icosahedron 

m 

12 

18 

21 

30 

Number of solutions 





complex 

72 

256 

12580 

- 

real (up to reflections) 

4 

18 

2 

- 

on S 

4 

8 

2 

32 

Running times 





HOM4PS-2.0 

19.8s 

3776s 

lOh 

- 

tracking 5 

0.06s 

0.8s 

0.07s 

8.6s 

Algo. 1^ 

0.02s 

0.09s 

0.01s 

1.5s 


7.3.2 Data of Table [T] 

Table gives for each problem the number m of equations of the system In group of lines 
“Number of solutions”, it first gives the total number of complex solutions of |J^so[ all found by 
HOM4PS-2.0. (line “complex”). The line “real (up to reflections)” gives the number of real 
solutions up to reflections. The line “on S” gives the number of real solutions lying on the path 
S to which belongs the sketch, that is obtained with our method. As remarked above, these 
solutions are different up to reflections. 

The group of lines “Running times” refers to times required to solve each problem with each 
approach. When using HOM4PS-2.0., most efforts are spent to follow paths leading to complex 
solutions, what explains the large running times in the line “HOM4PS-2.0.”. The line “tracking 
S” refers to the time required to track S in the space M"* x M, without using a RCP. The line 
“Algo. [2’ refers to the time required to compute solutions on S with on-the-fly change of RCP; 
it allows an important gain in term of computation cost. 

For the icosahedron problem, we do not give total number of solutions and running times for 
HOM4PS-2.0 since the solving process did not finish. 

7.3.3 Details on execution 

Table gives for each problem details about paths tracking without RCP (in the first columns) 
and with RCP (in the second columns). The row “time t in s” recalls the execution time in 
seconds. 

The row “smallest 6" gives the smallest prediction step used during the tracking process. It 
shows that tracking a path with a RCP requires to take smaller prediction steps than without 
a RCP. Together with the fact that Smax has to be smaller when using a RCP, it suggests that 
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Table 2: Details about path-tracking. For each problem, the first (resp. second) columns refers 
to the path tracking without (resp. with) RCP. 


m 

Octahedron 

12 

Disulfide 

18 

Dodecagon 

21 

Icosahedron 

30 

time t in s 

0.06 

0.02 

0.8 

0.09 

0.07 

0.01 

8.6 

1.5 

smallest 5 

0.1 

0.03 

0.1 

0.01 

0.05 

0.02 

0.05 

0.001 

nb. i of iterations 

157 

173 

1097 

694 

42 

46 

4074 

3377 

t/i in ms 

0.4 

0.1 

0.78 

0.13 

1.78 

0.3 

2.1 

0.43 

nb. of RCP changing 

- 

13 

- 

28 

- 

2 

- 

250 

average nb. of DP 

- 

1.3 

- 

1.98 

- 

4 

- 

3.8 

max. nb. of DP 

- 

2 

- 

3 

- 

4 

- 

6 


i?-reduced paths have higher curvature than the corresponding paths in the space of figures, and 
are more difficult to track. 

The lines “nb. i of iterations” and “t/i in ms” give the number of iterations of prediction- 
correction and the average time needed for each iteration. The latter information underlines the 
main advantage of using a RCP: each iteration of prediction-correction involves less computation 
than without RCP since the path is tracked in a space of smaller dimension. The line “time t in 
s” shows that this gain counterbalances the drawback mentioned above. 

The line “nb. of RCP changing” gives the number of times the RCP has been changed during 
the tracking process. This number can be large (see the case of the icosahedron problem): as 
mentioned in the penultimate paragraph of Sec. |3.4[ figures with a boundary configuration are 
the solutions of systems of to -|- 1 equations in to -I- 1 unknowns and can be in an exponential 
number. Hence the approach proposed in this paper could require, in the worst case, to change 
exponentially many times the RCP. The rows “average nb. of DP” and “max. nb. of DP” give 
respectively the average and the maximum number of driving parameters involved in the RCP 
and show that the number of driving parameters involved in successive RCP stays much lower 
than TO what keeps the method efficient. 


8 Conclusion 

Well-constrained point distance solving problems often have many solutions. The existing solvers 
that offer all the solutions are of limited practical interest because either the class of problems 
they solve is reduced or their complexity is exponential. But even if not all solutions are needed, 
several ones similar in shape to the sketch must be provided. 

An approach to fulfill this requirement is to use the sketch to define a real homotopy such that 
the homotopy path to which belongs the sketch is diffeomorphic to a circle and contains several 
solutions, that are similar to the sketch in the sense that they belongs to the same homotopy 
path. 

In this article we made this approach more efficient by reducing the dimension of the space 
where the homotopy path is tracked by using a symbolic geometric constructions program. The 
latter is modified on-the-fly in order to stay robust to critical geometric configurations it could 
induce. 

This original idea has been implemented to prove its soundness. In the examples discussed 
solutions are produced quicker when a construction program is used. Moreover the presented 
experiments show that our approach can provide several solutions to problems that are too large 
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to be solved with numerical solvers searching all the solutions such as homotopy. 

Notice finally that our method could be extended to more general geometric constraints 
such as angles, collinearities, coplanarities, and so on. When considering these constraints, 
homotopy paths are not necessarily diffeomorphic to circles but can converge to special geometric 
configurations that can be detected when using a construction program to stop the path tracking 
process. 
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A Numerical values of and 

Octahedron Values for parameters are picked up from [5]. 

Pi = (0,0,0), P2 = (1,0,0), P3 = (0.5, 1.3,0) 

P4 = (1.2, 0.5, 1.5), P5 = (0, 1.8,0.9), P6 = (-0.2,-0.7, 1.3) 

Disulfide molecule The reader is referred to m to get values of parameters. 

Pi = (0,0,0), P2 = (2.5, 5.3,0), P3 = (3.2, 5.8, 5.9), pi = (1, 8.7, 8.3), 

P5 = (4.2, 5.6, 7.9), P6 = (-2.2, 1.3, 6.3), P7 = (5.5, 2.5, 6.5), pg = (5,0,0). 


Dodecagone In the tables below, Oij denotes the value of the parameter of distance{pi,Pj). 


(21,2 — 3, <22,3 — 1.75, 

<27,8 = 1.45, as,9 — 1.35 

ai,4 = 4.4, 03,6 = 5.1, 

ai,7 = 4.45, as,9 = 6.65, 


<23,4 = 1.7, 
<29,10 = 1^ 
<25,8 — 3.9, 
<25,11 = 4.65. 


<24,5 = 2.05 
^10,11 = 1-4, 
07,10 = 3.05, 


<25,6 — 1-5, 
ail,12 = 1, 
09,12 = 3.35, 


06,7 — 1.85, 
012,1 — 0.6, 
011,2 = 4.4, 


Pi = (0,0), P2 = (2.9,0), 

P5 = (1.3, 4.6), P6 = (0.1, 5.7), 

P9 = (-3.0, 2.3), pio = (-2.4, 1.5), 


P3 = (3.4, 1.8), P4 = (2.8, 3.4), 
P7 = (-1.0,4.4), p8 = (-2.2,3.6), 
pii = (-1.2, 0.8), pi2 = (-0.6, 0.4) 


Icosahedron 


01,2 = 

2.0, 

02,4 = 

4, 

04,3 = 4:.5, 

03,6 = 

: 3.9, 

06,5 = 

4.45, 


05,8 — 3.8, 

08,9 = 

4.4, 

09,10 = 

^ 3.7, 

Oio,i2 = 4.35, 

012,11 

= 3.65, 

011,2 = 

- 4.3, 


Oil,7 — 2.5, 

04,7 = 

2.6, 

06,7 = 

2.7, 

as,7 — 2.8, 

Ol0,7 - 

- 2.9, 

Oil,4 = 

= 3.0, 


Oil,8 — 2.95, 

Ol0,8 = 

= 3.05, 03,6 = 

2.9, 

06,4 =3.1, 

02,3 = 

: 3.0, 

03,5 — 

2.9, 


05,9 — 3.1, 

09,12 = 

= 2.8, 

Ol2,2 = 

= 3.2, 

Ol,9 = 2.1, 

Oi,5 = 

: 2.2, 

Ol,3 = 

2.3, 


01,12 — 2.4. 

Pi 

= (0 

.0,0), 

P2 = 

= (1.0,0), 

P3 ^ 

(0.5,1,0), 


P4 ^ 

(IL 

1.2, 1.3), 

P5 

= (- 

-0.5,0.5,1), 

P6 = 

= (1.7,2.3.1.1), 

P7 ^ 

(3,3,3), 


P8 ^ 

(0.8 

, 1.6,2.2), 

P9 

= (- 

-0,0,2), 

PlO 

= ( 1 . 1 , 1 . 3 , 3 . 2 ), 

Pll = 

= ( 1 . 8 , 1 . 2 , 

2 . 1 ), 

P12 = 

= (0.1 

6,0.2, 1.3) 


Inria 




RESEARCH CENTRE 
NANCY-GRAND EST 

615 rue du Jardin Botanique 
CS20101 

54603 Villers-les-Nancy Cedex 


Publisher 

Inria 

Domaine de Voluceau - Rocquencourt 
BP 105 - 78153 Le Chesnay Cedex 
inria.fr 

ISSN 0249-6399 



